################################################################################
#                                                                              #
# TITLE        : BREAKING SILENCE: HOW INTIMATE PARTNER VIOLENCE AND           #         
#                REPORTING SHAPE LATER LIFE OUTCOMES                           #        
# CODE EDITOR  : Harrison Chang                                                #
# LAST MODIFIED: 2025/07/18                                                    #
# PURPOSE      : This code explores reporting behavior                         #
#                                                                              #
#                                                                              #
#                                                                              #
#                                                                              #
################################################################################


rm(list=ls());gc()

library(data.table)
library(stringr)
library(tidyverse)
library(dplyr)
library(fixest)
library(nnet)
library(broom)
library(knitr)
library(kableExtra)
library(gtsummary)
library(ggeffects)
library(marginaleffects)
library(car)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
####         Table B.6: Multinomial Logit Model      ####
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

#################################################
#            Data Reading and Cleaning          # 
#################################################

load("E:/H112057/H112057_Harrison/DV/rdata/IPVf_once_report_all.Rdata")

unique(TMP, by = "ID_no")[, .N, .(event_year)][order(event_year)]
unique(TMP, by = "ID_no")[, .N, .(infounit)][order(infounit)]

#Initially, it is 1(D)=Police; 2(C)=Hospital; 3(I)=Helpline; 4(other)=other 
#But now I want to switch D and C
TMP[infounit==2, temp:=1]
TMP[infounit==1, temp:=2]
TMP[temp==1    ,infounit:=1]
TMP[temp==2    ,infounit:=2]
unique(TMP, by = "ID_no")[, .N, .(infounit)][order(infounit)]
#So now it's 1(C)=Hospital; 2(D)=Police; 3(I)=Helpline; 4(other)=other 


#First generate relative year, then only keep those one year beforehand
TMP[, relative_year := calendarYr - event_year]
TMP[, .N, .(relative_year)][order(relative_year)]
cross_section   <- TMP[relative_year==-1][!duplicated(ID_no)]

cross_section[, vic_age := event_year - birth_year]
cross_section[, .N, .(vic_age)][order(vic_age)]

#age group
cross_section[vic_age %in% c(25:35), age_group := 1]
cross_section[vic_age %in% c(36:45), age_group := 2]
cross_section[vic_age %in% c(46:55), age_group := 3]
cross_section[, .N, .(age_group)][order(age_group)]


#edu group
cross_section[, .N, .(max_edu_year)][order(max_edu_year)]
cross_section[max_edu_year <= 9, edu_group := 1]
cross_section[max_edu_year ==12, edu_group := 2]
cross_section[max_edu_year >=14, edu_group := 3]
cross_section[, .N, .(edu_group)][order(edu_group)]

#child number group
cross_section[, .N, .(num_children_under_12)][order(num_children_under_12)]

#victim income group
cross_section[, .N, .(TRUE_AMT)][order(TRUE_AMT)]
summary(cross_section$TRUE_AMT)
cross_section[floor(TRUE_AMT) ==0                   , income_group := 1]
cross_section[floor(TRUE_AMT) %in% c(1:264759)      , income_group := 2]
cross_section[floor(TRUE_AMT) %in% c(264760:305671) , income_group := 3]
cross_section[floor(TRUE_AMT) %in% c(305672:2295324), income_group := 4]
cross_section[, .N, .(income_group)][order(income_group)]

cross_section[infounit==1, channel := "Hospital"]
cross_section[infounit==2, channel := "Police"]
cross_section[infounit==3, channel := "Helpline"]
cross_section[infounit==4, channel := "Other"]

cross_section$channel <- cross_section$channel %>% as.factor()
cross_section$channel <- relevel(cross_section$channel, ref = "Hospital")




#################################################
#                Panel A - Vic_Age              # 
#################################################

mymodel    <- cross_section %>% copy()

test       <- multinom(channel ~ vic_age +
                                 edu_group + 
                                 income_group +
                                 num_children_under_12, data = mymodel)

###First show the change in log odds
summary(test)
tidy(test)



mfx <- slopes(test, variables = "vic_age") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group <- mfx %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat <- mean_se_by_group$mean_marginal / mean_se_by_group$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value <- 2 * (1 - pnorm(abs(t_stat)))  # Two-tailed test
# View the results
export <- cbind(mean_se_by_group,p_value)
export



#################################################
#           Panel A - Education Group           # 
#################################################


mfx2 <- slopes(test, variables = "edu_group") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group2 <- mfx2 %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat2 <- mean_se_by_group2$mean_marginal / mean_se_by_group2$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value2 <- 2 * (1 - pnorm(abs(t_stat2)))  # Two-tailed test
# View the results
export2 <- cbind(mean_se_by_group2,p_value2)
export2






#################################################
#          Panel A - NO. Child Under 12         # 
#################################################


mfx3 <- slopes(test, variables = "num_children_under_12") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group3 <- mfx3 %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat3 <- mean_se_by_group3$mean_marginal / mean_se_by_group3$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value3 <- 2 * (1 - pnorm(abs(t_stat3)))  # Two-tailed test
# View the results
export3 <- cbind(mean_se_by_group3,p_value3)
export3







#################################################
#              Panel A - Income Group           # 
#################################################


mfx4 <- slopes(test, variables = "income_group") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group4 <- mfx4 %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat4 <- mean_se_by_group4$mean_marginal / mean_se_by_group4$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value4 <- 2 * (1 - pnorm(abs(t_stat4)))  # Two-tailed test
# View the results
export4 <- cbind(mean_se_by_group4,p_value4)
export4







#################################################
#     Panel B - Perpetrator Education Group     # 
#################################################


cross_section <- cross_section[is.na(perp_edu)==0]

#perp edu group
cross_section[, .N, .(perp_edu)][order(perp_edu)]
cross_section[perp_edu <= 9 , perp_edu_group := 1]
cross_section[perp_edu ==12 , perp_edu_group := 2]
cross_section[perp_edu >=14 , perp_edu_group := 3]
cross_section[, .N, .(perp_edu_group)][order(perp_edu_group)]


#perp income group
cross_section[, .N, .(perp_TRUE_AMT)][order(perp_TRUE_AMT)]
summary(cross_section$perp_TRUE_AMT)
cross_section[floor(perp_TRUE_AMT) ==0                   , perp_income_group := 1]
cross_section[floor(perp_TRUE_AMT) %in% c(1:271101)      , perp_income_group := 2]
cross_section[floor(perp_TRUE_AMT) %in% c(271102:400866) , perp_income_group := 3]
cross_section[floor(perp_TRUE_AMT) %in% c(400867:2295324), perp_income_group := 4]
cross_section[, .N, .(perp_income_group)][order(perp_income_group)]

mymodel2      <- cross_section %>% copy()


test_perp     <- multinom(channel ~ vic_age +
                                    edu_group + 
                                    income_group +
                                    num_children_under_12 +
                                    perp_edu_group +
                                    perp_income_group, data = mymodel2)



mfx5 <- slopes(test_perp, variables = "perp_edu_group") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group5 <- mfx5 %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat5 <- mean_se_by_group5$mean_marginal / mean_se_by_group5$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value5 <- 2 * (1 - pnorm(abs(t_stat5)))  # Two-tailed test
# View the results
export5 <- cbind(mean_se_by_group5,p_value5)
export5



#################################################
#       Panel B - Perpetrator Income Group      # 
#################################################


mfx6 <- slopes(test_perp, variables = "perp_income_group") %>% data.table()
# Calculate mean and standard error of AME for each group
mean_se_by_group6 <- mfx6 %>%
  group_by(group) %>%
  summarise(
    mean_marginal = mean(estimate),
    se_marginal = sd(std.error) / sqrt(n())  # Standard error: SD / sqrt(n)
  )
# Calculate the t-statistic manually
t_stat6 <- mean_se_by_group6$mean_marginal / mean_se_by_group6$se_marginal # AME divided by standard error
# Compute the p-value from the t-statistic
p_value6 <- 2 * (1 - pnorm(abs(t_stat6)))  # Two-tailed test
# View the results
export6 <- cbind(mean_se_by_group6,p_value6)
export6










# # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
####           Table B.8: IPV Duration Analysis      ####
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # 


#################################################
#           Reporting Behavior: Setting         #
#################################################

load("E:/H112057/H112057_Harrison/DV/rdata/IPVf_once_violence_all.Rdata")

unique(TMP, by = "ID_no")[, .N, .(event_year)][order(event_year)]
unique(TMP, by = "ID_no")[, .N, .(duration_year)][order(duration_year)]

#First generate relative year, then only keep those one year beforehand
TMP[, relative_year := calendarYr - event_year]
TMP[, .N, .(relative_year)][order(relative_year)]
cross_section   <- TMP[relative_year==-1][!duplicated(ID_no)]

#sum stats
unique(cross_section, by = "ID_no")[, .N, .(duration_year)][order(duration_year)]
cross_section[, .N, .(event_year,duration_year)][order(event_year,duration_year)]

#We generate the first set of covariate - victimized age
cross_section[, vic_age := event_year - birth_year]
cross_section[, .N, .(vic_age)][order(vic_age)]

#Then the second set - victim income dummy
cross_section[, TRUE_AMT := floor(TRUE_AMT)]
summary(cross_section$TRUE_AMT)

cross_section[, income_group4 := ifelse(TRUE_AMT==0,1,0)]
cross_section[, income_group3 := ifelse(TRUE_AMT %in% c(1:258000),1,0)]
cross_section[, income_group2 := ifelse(TRUE_AMT %in% c(258001:302950),1,0)]
cross_section[, income_group1 := ifelse(TRUE_AMT %in% c(302951:2295324),1,0)]
cross_section[, .N, .(income_group1)]
cross_section[, .N, .(income_group2)]
cross_section[, .N, .(income_group3)]
cross_section[, .N, .(income_group4)]


#Then the second set - victim income dummy
cross_section[, perp_TRUE_AMT := floor(perp_TRUE_AMT)]
summary(cross_section$perp_TRUE_AMT)

cross_section[, perp_income_group4 := ifelse(perp_TRUE_AMT==0,1,0)]
cross_section[, perp_income_group3 := ifelse(perp_TRUE_AMT %in% c(1:270955),1,0)]
cross_section[, perp_income_group2 := ifelse(perp_TRUE_AMT %in% c(270956:382133),1,0)]
cross_section[, perp_income_group1 := ifelse(perp_TRUE_AMT %in% c(382134:2295324),1,0)]
cross_section[, .N, .(perp_income_group1)]
cross_section[, .N, .(perp_income_group2)]
cross_section[, .N, .(perp_income_group3)]
cross_section[, .N, .(perp_income_group4)]



#Approximate victim ability
cross_section[, max_pre_income := floor(max_pre_income)]
summary(cross_section$max_pre_income)

cross_section[, max_income_group4 := ifelse(max_pre_income %in% c(0:232353),1,0)]
cross_section[, max_income_group3 := ifelse(max_pre_income %in% c(232354:276195),1,0)]
cross_section[, max_income_group2 := ifelse(max_pre_income %in% c(276196:387550),1,0)]
cross_section[, max_income_group1 := ifelse(max_pre_income %in% c(387551:2295324),1,0)]
cross_section[, .N, .(max_income_group1)]
cross_section[, .N, .(max_income_group2)]
cross_section[, .N, .(max_income_group3)]
cross_section[, .N, .(max_income_group4)]


#Approximate perpetrator ability
cross_section[, perp_max_pre_income := floor(perp_max_pre_income)]
summary(cross_section$perp_max_pre_income)

cross_section[, max_perp_income_group4 := ifelse(perp_max_pre_income %in% c(0:247955),1,0)]
cross_section[, max_perp_income_group3 := ifelse(perp_max_pre_income %in% c(247956:289064),1,0)]
cross_section[, max_perp_income_group2 := ifelse(perp_max_pre_income %in% c(289065:490135),1,0)]
cross_section[, max_perp_income_group1 := ifelse(perp_max_pre_income %in% c(490136:2295324),1,0)]
cross_section[, max_perp_income_group5 := ifelse(is.na(perp_max_pre_income)==1,1,0)]
cross_section[, .N, .(max_perp_income_group1)]
cross_section[, .N, .(max_perp_income_group2)]
cross_section[, .N, .(max_perp_income_group3)]
cross_section[, .N, .(max_perp_income_group4)]
cross_section[, .N, .(max_perp_income_group5)]


#perpetrator education
cross_section[, .N, .(perp_edu)][order(perp_edu)]

#report_source
cross_section[, .N, .(INFOUNIT)]
cross_section[, source4 := ifelse(INFOUNIT=="C",1,0)]
cross_section[, source2 := ifelse(INFOUNIT=="D",1,0)]
cross_section[, source3 := ifelse(INFOUNIT=="I",1,0)]
cross_section[, source1 := ifelse(INFOUNIT!="C"&INFOUNIT!="D"&INFOUNIT!="I",1,0)]
cross_section[, .N, .(source1)]
cross_section[, .N, .(source2)]
cross_section[, .N, .(source3)]
cross_section[, .N, .(source4)]

#occupation
cross_section[, private_worker := ifelse(occupation=="B",1,0)]


#################################################
#         Reporting Behavior: Regression        #
#################################################

model0          <- feols(duration_year ~ vic_age +
                                         max_edu_year + 
                                         num_children_under_12, data = cross_section)


model1          <- feols(duration_year ~ vic_age +
                                         max_edu_year + 
                                         num_children_under_12 | event_year , data = cross_section)

model2          <- feols(duration_year ~ vic_age +
                                         max_edu_year + 
                                         num_children_under_12 +
                                         income_group1 + 
                                         income_group2 +
                                         income_group3 +
                                         income_group4 +
                                         year_full_dummy | event_year, data = cross_section)

model3          <- feols(duration_year ~ vic_age +
                                         max_edu_year + 
                                         num_children_under_12 +
                                         income_group1 + 
                                         income_group2 +
                                         income_group3 +
                                         income_group4 +
                                         year_full_dummy +
                                         source1 +
                                         source2 +
                                         source3 +
                                         source4 | event_year, data = cross_section)

cross_section   <- cross_section[max_perp_income_group5!=1]

model4          <- feols(duration_year ~ vic_age +
                                         max_edu_year + 
                                         num_children_under_12 +
                                         income_group1 + 
                                         income_group2 +
                                         income_group3 +
                                         income_group4 +
                                         year_full_dummy +
                                         source1 +
                                         source2 +
                                         source3 +
                                         source4 + 
                                         perp_edu +
                                         perp_income_group1 +
                                         perp_income_group2 +
                                         perp_income_group3 +
                                         perp_income_group4 | event_year, data = cross_section)




etable(model1,
       model2,
       model3,
       model4,
       headers = c("(1)","(2)","(3)","(4)"))













